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THE STELLAR POPULATION STRUCTURE OE THE GALACTIC DISK 
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Matthew Shetrone^, and Timothy C. Beers® 

ABSTRACT 

The spatial structure of stellar populations with different chemical abundances in the Milky Way 
contains a wealth of information on Galactic evolution over cosmic time. We use data on 14,699 
red-clump stars from the APOGEE survey, covering 4kpc ^ R ^ 15kpc, to determine the structure 
of mono-abundance populations (MAPs)—stars in narrow bins in [o/Ee] and [Ee/H]—accounting for 
the complex effects of the APOGEE selection function and the spatially-variable dust obscuration. 
We determine that all MAPs with enhanced [a/Ee] are centrally concentrated and are well-described 
as exponentials with a scale length of 2.2 ± 0.2 kpc over the whole radial range of the disk. We 
discover that the surface-density profiles of low-[(T/Ee] MAPs are complex: they do not monotonically 
decrease outwards, but rather display a peak radius ranging from ^ 5kpc to ^ 13kpc at low [Ee/H]. 
The extensive radial coverage of the data allows us to measure radial trends in the thickness of each 
MAP. While high-[(T/Ee] MAPs have constant scale heights, low-[(T/Ee] MAPs flare. We confirm, 
now with high-precision abundances, previous results that each MAP contains only a single vertical 
scale height and that low-[Ee/H], low-[(T/Ee] and high-[Ee/H], high-[(T/Ee] MAPs have intermediate 
{hz ~ 300-600 pc) scale heights that smoothly bridge the traditional thin- and thick-disk divide. That 
the high-[(T/Ee], thick disk components do not flare is strong evidence against their thickness being 
caused by radial migration. The correspondence between the radial structure and chemical-enrichment 
age of stellar populations is clear confirmation of the inside-out growth of galactic disks. The details 
of these relations will constrain the variety of physical conditions under which stars form throughout 
the MW disk. 

Subject headings: Galaxy: abundances — Galaxy: disk — Galaxy: evolution — Galaxy: formation 
— Galaxy: fundamental parameters — Galaxy: structure 


1. INTRODUCTION 

Understanding the growth and evolution of galactic 
disks is a central problem in galaxy formation. Em¬ 
pirical studies in this area employ two complementary 
approaches to investigate changes in the structure of 
galactic disks over cosmic time. One of them, the “look- 
back approach”, is based on extensive samples of galax¬ 
ies that span redshifts from the time at which disks 
form {z ~ 2.5) to today {z ~ 0), matching correspond¬ 
ing samples across epochs (e.g., van Dokkum et al. 2013; 
Wisnioski et al. 2015). Such studies have been success¬ 
ful in establishing the evolution of the galaxy population 
properties, but they are limited to characterizing indi¬ 
vidual galaxies by quantities integrated over their stellar 
populations, e.g., by their size, shape, or overall veloc¬ 
ity dispersion. The second approach. Galactic archae¬ 
ology, aims at obtaining detailed observations of local 
galaxies, either through resolving their stellar popula- 
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tions into individual stars (e.g., Dalcanton et al. 2012) or 
using integral-field spectroscopy, trying to directly recon¬ 
struct their individual formation history. The Milky Way 
(MW) is perhaps the best case for Galactic archaeology, 
because it has very typical gross properties (e.g., mass) 
and because we can measure the detailed properties— 
six-dimensional phase-space distribution, age, elemen¬ 
tal abundances—of large numbers of stars (Rix & Bovy 
2013). 

The growth of the MW disk over time is encoded in 
the orbital distribution of stars and their ages and el¬ 
emental abundances. The radial distribution of stars 
of different ages and abundances constrains in partic¬ 
ular the epochs and physical conditions under which 
different parts of the MW formed stars. However, the 
large-scale radial structure of the MW disk has been dif¬ 
ficult to determine for any well-defined stellar tracer, 
primarily due to severe dust extinction in the mid¬ 
plane. Investigations have therefore been limited to 
measurements of the overall radial profile, which can 
be characterized as an approximately exponential disk 
(e.g., Kent et al. 1991; Benjamin et al. 2005; Juric et al. 
2008; Bovy & Rix 2013). Radial population changes 
in the Galactic disk have traditionally been charac¬ 
terized by a metallicity gradient, taking at each ra¬ 
dius the mean over the complex abundance distribu¬ 
tion (e.g., Audouze & Tinsley 1976; Chen et al. 2003; 
Anders et al. 2014; Hayden et al. 2014). Neither of these 
measurements provides a direct empirical constraint on 
the growth of the MW disk over time, nor have they 
proven to be very constraining for simulations of the for¬ 
mation of galactic disks. 
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Fig. 1.— Distribution of the APOGEE fields that contain the statistical APOGEE-RC sample on the sky, overlayed on the extinction 
map from Green et al. (2015) at 5kpc. Fields at |6| < 8° are displayed in white to enhance their contrast. The upper limit of Ah = 0.8 mag 
is the highest extinction for which the RG can be seen to 5kpc in APOGEE medium-deep fields (H < 12.2). The statistical APOGEE-RG 
has excellent coverage of the large portion of the Galactic plane that can be seen from the Northern hemisphere (the large white region is 
unobserved by both APOGEE and Pan-STARRS). 


Bovy et al. (2012c) suggested a different way to look 
at the complex correlations between spatial structure 
and abundances of the Galactic disk population, by sep¬ 
arately determining disk structure of so-called mono¬ 
abundance populations (MAPs). These are subsets of 
stars selected to have very similar abundances, such as 
[Fe/H] and and average [a/Fe]; for each star these abun¬ 
dances constitute life-long tags that remain invariant 
even if their orbits change grossly after birth. Beyond 
this simple fact, a decomposition of the disk in terms 
of MAPs is essentially empirical and does not assume 
or imply a common formation site or epoch; the detailed 
relation between MAPs and the underlying chemical evo¬ 
lution of the MW requires models and measurements of 
the age distributions within MAPs. Bovy et al. (2012c) 
were the first to attempt dissecting the radial stellar 
population structure of the MW disk, by measuring the 
radial profiles of MAPs over a few kpc near the Sun. 
These measurements revealed a complex abundance- 
dependence of the radial disk structure—both on [Fe/H] 
and an average [a/Fe]—with older, [a/Fe]-enhanced 
populations having centrally-concentrated profiles and 
younger, solar-[a/Fe] populations having more extended 
profiles reaching the outer disk. However, the limited 
radial coverage and high-latitude survey geometry of the 
SDSS/SEGUE sample (Yanny et al. 2009) impeded more 
detailed measurements of the radial surface density pro¬ 
file beyond characterizing the local slope (in radius) by 


an exponential-disk scale length. 

The vertical profile of the MW, near the Sun and 
at other radii, is also of great interest for understand¬ 
ing the MW, although its interpretation in terms of 
stellar disk growth and evolution is still unsettled. In 
the solar neighborhood, the overall vertical stellar den¬ 
sity profile (that is, of all stars) can be accurately 
characterized as a sum of two exponential components 
(Yoshii 1982; Gilmore & Reid 1983; Juric et al. 2008; 
Rix & Bovy 2013), with stars in the thicker of the two 
components being older, more metal-poor, and [a/Ee]- 
enhanced with respect to stars in the thinner components 
(Euhrmann 1998; Prochaska et al. 2000; Bensby et al. 
2005; Reddy et al. 2006; Haywood et al. 2013). Dissect¬ 
ing the disk into MAPs, Bovy et al. (2012a,b,c) demon¬ 
strated that this broader picture is the consequence of 
an underlying disk structure composed of a continuum 
of disks, with a smooth correlation between abundance, 
scale height, and velocity dispersion. While the vertical 
profile near the Sun is composed of a smooth contin¬ 
uum of disk thicknesses, it has become clear from unbi¬ 
ased observations of element abundances of disk stars 
that their distribution in ([Ee/H], [a/Ee]) is distinctly 
bimodal, composed of two sequences—high- and low- 
[a/Ee]—that are disjoint at low-[Ee/H] and that merge 
near solar [Ee/H] (Adibekyan et al. 2012; Nidever et al. 
2014). The high-[(a/Ee] sequence becomes more promi¬ 
nent in the inner Galaxy, while the low-[(a/Ee] and in par- 






3 




Fig. 2.— Peak of the absolute-magnitude PDF p(Mh\[J — Ks]o) 
as a function of color (J — Ks)o and metallicity Z for PARSEC 
isochrones in the RC-star region (defined using Equations [2] and 
[3] of Bovy et al. 2014). The white dashed lines represent the region 
specified by the cuts in Equations (6) and (7) of Bovy et al. (2014) 
that delineate the regions over which the distribution of absolute 
Ks magnitudes is narrow; these also select regions where Mh is 
narrowly distributed. The peak of the magnitude PDF does not 
strongly depend on color or metallicity over the region where the 
PDF is narrow. These Mh{[J — Ks]o,Z) are combined with an 
overall offset of AM^ = 0.08 to determine the distances for RC 
stars used in this paper. 

ticular its metal-poor end dominate in the outer Galaxy 
(Bensby et al. 2011; Nidever et al. 2014; Hayden et al. 
2015), in agreement with the scale lengths measured from 
the SEGUE data (Bovy et al. 2012c). But a more de¬ 
tailed and quantitative characterization of the radial be¬ 
havior of the stellar populations does not yet exist. Ex¬ 
actly how the continuity in local disk thicknesses and the 
disjoint high- and low-[a/Ee] sequences in abundances 
can be reconciled in a consistent formation scenario is 
currently a mystery. The complex behavior of the ra¬ 
dial scale lengths with abundance found by Bovy et al. 
(2012c), as opposed to the smooth variation of the verti¬ 
cal scale height, indicates that a thorough exploration of 
the radial profile of MAPs might hold the key to resolve 
this tension. 

While the MW disk’s overall vertical stellar den¬ 
sity profile was long thought to be the result of 
some type of heating, whether from satellites (e.g., 
Quinn et al. 1993; Abadi et al. 2003; Brook et al. 2004) 
or through radial migration (e.g., Sellwood & Binney 
2002; Schonrich & Binney 2009; Minchev & Eamaey 
2010; Loebman et al. 2011), the possibility that (at least 
some part of) the stellar disk formed “upside-down” 
with its current thickness has recently gained favor 
(Bournaud et al. 2009; Stinson et al. 2013; Bird et al. 

2013) . The smooth continuum of disk thicknesses found 
by Bovy et al. (2012a,c) disfavors the satellite accretion 
and few-satellites-heating scenarios, because those would 
create a discrete thick-disk component (e.g., Martig et al. 

2014) . However, observational evidence to prefer either 
radial migration or upside-down formation is scant to 
date. 

Radial migration of stars from well inside the MW 
has been proposed as a mechanism to create the 
thicker, [a/Ee]-enhanced stellar disk components in 


the solar neighborhood. Radial migration approxi¬ 
mately conserves the vertical action (Solway et al. 2012; 
Vera-Giro & D’Onghia 2015). Eor outwardly migrating 
stars, the decreased gravitational pull in the outer disk 
leads to larger vertical excursions that for ensembles of 
stars should lead to flaring. Whether radial migration 
causes the whole stellar disk or its components to flare 
depends on which stars participate in the migration pro¬ 
cess. If all stars migrate similarly, then all co-eval pop¬ 
ulations should flare, but if only special subclasses are 
affected (e.g., stars with low vertical excursions), then 
the overall effect of migration on the vertical structure 
might be small. While the current suite of simulations by 
no means exhaustively covers the various migration sce¬ 
narios and histories that may have occurred in the MW, 
current simulations do indicate that kinematic biases in 
the population of migrators are significant enough that 
little or no effect on the vertical structure of the disk 
should be expected (Minchev et al. 2012; Roskar et al. 
2013; Vera-Ciro et al. 2014). Significant radial mixing, 
therefore, may well happen without leading to flaring 
while still having a large effect on, e.g., the present-day 
abundance distribution. However, for migration to be 
the cause of the thickness of the thicker-disk compo¬ 
nents in the MW, flaring of sub-populations must occur 
and determining whether individual MAPs flare verti¬ 
cally towards larger radii can discriminate among forma¬ 
tion mechanisms for the thick-disk components. 

That radial migration should and does play a large 
role in the evolution of the low-[a/Ee] populations is be¬ 
coming increasingly clear. Radial migration provides a 
natural way to explain the observed lack of correlation 
between age and metallicity in the solar neighborhood 
(Edvardsson et al. 1993; Nordstrom et al. 2004) in the 
presence of a strong correlation between age and [a/Ee] 
(Haywood et al. 2013). The change in skew of the metal¬ 
licity distribution from negative in the inner MW to posi¬ 
tive in the outer MW (Hayden et al. 2015) can arise from 
migration. As for the thicker components, if this migra¬ 
tion process is important for all low-[(a/Ee] stars (not 
just a biased subset), it makes the generic prediction of 
an outward flaring disk. 

In this paper we set out to comprehensively de¬ 
termine the stellar population structure of the Milky 
Way’s stellar disk. We do so by determining the radial 
and vertical profile of MAPs, using SDSS-HI/APOGEE 
(Eisenstein et al. 2011; Majewski et al. 2015) data. The 
APOGEE data are complementary to the SEGUE data, 
because they provide much better radial coverage by us¬ 
ing near-infrared observations of a large number of gi¬ 
ants close to the Galactic midplane: the sample of red- 
clump giants used here spans the radial range 4kpc < 
R < 15kpc. This radial coverage allows us to determine 
the radial profile of MAPs in detail: the data clearly re¬ 
quire to go beyond the simple radial-exponential profile 
assumed by Bovy et al. (2012c) and all other previous 
analyses. The wide radial range combined with good 
vertical coverage (0 < |Z| < 2kpc) also makes it possible 
to measure any flaring of MAPs, providing an essential 
test of the predictions of radial migration. Einally, the 
high-resolution spectroscopic observations let us define 
MAPs with high-precision ([Ee/H], [a/Ee]) abundances, 
such that MAPs have no substantial contamination from 
neighboring MAPs. These better APOGEE data end up 
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Fig. 3.— Fractional difference between the RC distances calcu¬ 
lated based on the if-band and Kg-hand luminosities as a function 
of the latter distance. The conditional distribution of distance 
differences is displayed for all 14,699 stars in the APOGEE-RC 
subsample used in this paper. The black contours contain 68 % 
and 95 % of the distribution and the white dashed line shows a 
lowess trendline. The overall median difference is only 2.3%, con¬ 
sistent with the estimated '^2% systematic uncertainty in the RC 
distance scale (Bovy et al. 2014). 


confirming the smooth, continuous distribution in scale 
height from the thinnest to the thickest disk components, 
as found by Bovy et al. (2012a). 

The outline of this paper is as follows. In § 2 we 
present and discuss important aspects of the data from 
the APOGEE-RC catalog that this analysis is based on, 
paying particular attention to the abundance measure¬ 
ments in § 2.2. We describe the method for fitting the 
density profiles of stellar subsamples in APOGEE in § 3; 
this method is presented in more detail in Bovy et al. 
(2016). We apply the density-fitting methodology to 
broad abundance-selected subsamples in § 4 to explore 
the density profiles that represent the data well. In § 5 we 
then determine the density profiles of MAPs. We com¬ 
pare our results to previous work and discuss implications 
for our understanding of the formation and evolution of 
the MW disk in § 6. Einally, we present our conclusions 
in § 7. Readers not interested in the details of the data 
and fitting methodology are encouraged to skip to § 4. 
Throughout this paper, we assume that the Sun’s dis¬ 
placement from the mid-plane is 25 pc toward the North 
Galactic Pole (Chen et al. 2001; Juric et al. 2008), and 
that the Sun is located at 8 kpc from the Galactic center. 


2. DATA 

2.1. APOGEE and the APOGEE-RC catalog 

The SDSS-III/APOGEE (Majewski et al. 2015) is a 
high-resolution {R^ 22, 500) spectroscopic survey in the 
near-infrared (NIR; i7-band; 1.51 to 1.70 jim). The sur¬ 
vey employs a 300-fiber spectrograph (Wilson et al. 2010, 


J. Wilson et al. 2016, in preparation) to obtain signal-to- 
noise ratio of 100 per half-resolution element (~ 141 per 
resolution element) spectra for large numbers of stars in 
the Milky Way, observed during bright time on the 2.5- 
meter Sloan Eoundation telescope (Gunn et al. 2006). 
The majority of APOGEE targets are selected from the 
2MASS point-source catalog (Skrutskie et al. 2006), us¬ 
ing a dereddened {J — Kg) o > 0.5 color cut in up to three 
magnitude bins in H (not corrected for extinction), with 
reddening corrections determined using the Rayleigh 
Jeans Color Excess method (RJCE; Majewski et al. 
2011) applied to 2MASS and mid-IR data from Spitzer- 
IRAC GLIMPSE-I, -II, and -3D (Churchwell et al. 2009) 
when available and from WISE (Wright et al. 2010) 
otherwise. In this paper, we employ various three- 
dimensional extinction maps to take the effects of ex¬ 
tinction on the observed iJ-band magnitude into ac¬ 
count, but we always use the RJCE dereddenings for 
the J — Kg color. A direct comparison between the 3D 
extinction maps and the RJCE extinctions exhibits good 
agreement between the two (see Bovy et al. 2016). Eull 
details on the APOGEE target selection can be found 
in Zasowski et al. (2013). We use data from the SDSS- 
III Data Release 12 (DR12; Alam et al. 2015). The 
data processing (Nidever et al. 2015), stellar-parameters 
and chemical-abundances analysis (Shetrone et al. 2015; 
Zamora et al. 2015; Garcia Perez et al. 2015, and the 
DR12 data analysis and calibration (Holtzman et al. 
2015) are performed with automated SDSS software; 
we make use of standard data products as contained in 
DR12. 

Specifically, we base our analysis on red-clump (RC) 
star data from the DR12 APOGEE-RC catalog. This 
catalog consists of a pure sample of RC stars selected 
from the APOGEE catalog using the method described 
in detail by Bovy et al. (2014), but applied to the DR12 
data (see Alam et al. 2015). RC stars are identified 
in the superset of all APOGEE data using a combi¬ 
nation of cuts in the surface-gravity (log^) / effective- 
temperature (Teff) plane and the {J — Kg) o / metallic- 
ity ([Ee/H]) plane, chosen such as to select those RC 
stars for which precise luminosity distances can be deter¬ 
mined. The catalog has an estimated purity of ~ 95 %; 
the distances are precise to 5 % and unbiased to 2 %. 
In order to account for biases due to the APOGEE tar¬ 
get selection, we only use the subset of RC data that 

(a) was selected as part of the “main” APOGEE sam¬ 
ple (defined using the [J — Kg]^ > 0.5 color cut), and 

(b) are part of a line of sight and magnitude bin com¬ 
bination for which the planned APOGEE observations 
were complete in DR12, because these are the stars 
for which we can reconstruct the sample selection func¬ 
tion. In our density-fitting methodology (see § 3 below) 
we correct for the effects of interstellar extinction us¬ 
ing the 3D extinction maps of Marshall et al. (2006) in 
the inner-disk plane and Green et al. (2015) elsewhere 
(using AnfAx^ = 1.48 and Ah/E{B — V) = 0.46; 
Schlafiy & Einkbeiner 2011; Yuan et al. 2013); the lat¬ 
ter is based on Pan-STARRS and 2MASS data and does 
not cover the entire APOGEE footprint. We there¬ 
fore remove 10 APOGEE fields that lie outside of the 
Green et al. (2015) footprint. These are the fields cen¬ 
tered on (/,6) = (240,-18), (5.5,-14.2), (5.2,-12.2), 
(1,-4), (0,-5), (0,-2), (358,0), (358.6,1.4) (difficult 
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Fig. 4.— Directly apparent signatures of [a/Fe] variations among APOGEE-RC stars, split by individual alpha-elements. This figure 
illustrates spectral differences in wavelength regions with strong features due to O, Mg, Si, S, Ca, and Ti for stars over a narrow range 
in [Fe/H] (—0.45 < [Fe/H] < —0.35). Each panel displays residuals from a mean spectrum after interpolating each spectral pixel to a 
common Tgff, log^f, and [Fe/H] using quadratic interpolation. Residuals are smoothed by only including the part of the spectrum contained 
in the sub-space spanned by the eigenvectors of the eight largest PC A components of the residuals of all 490 stars in this [Fe/H] range; 
this sub-space contains all of the variance above the measurement noise (we only include spectra with S/N > 200 in this figure). A 
stack of twelve random residuals each in five bins in A[Q;/Fe] = [([O + Mg + Si + S + Ca]/5)/Fe] = 0.05—which does not include Ti (see 
text)—ranging from 0.00 (dark blue) to 0.20 (dark red) is displayed. The x axis only covers small parts of the full APOGEE spectral range 
in each panel and is interrupted in most panels to focus on features due to a specific element; the wavelength of the reddest tickmark in 
each section is indicated and the tickmark spacing is 2 A everywhere. The location of the spectral features for each element are indicated 
in each panel. Note that the scale for the Ga panel is different due to the weakness of the Ga feature. This figure demonstrates the high 
precision in [a/Fe] and the high level of consistency between the abundances of different a elements obtained from the if-band APOGEE 
spectra: (^[a/Fe] ^ 0.02 dex. 


to access from the Northern Hemisphere) and (/,6) = 
(120,30), (123,22.4). This only removes 125 stars. We 
further remove 13 stars with distance moduli based on 
their i^-band luminosity (see below) < 8.49, because 
they could not be in our sample if we model the RC as 
a standard candle with Mh = —1.49 (our standard as¬ 
sumption below). This statistical sample contains 14,699 
RC stars. This sample spans a range 4500 K < Teff < 
5200K and 2.25 ^ log^f < 3 (95% ranges). 

Figure 1 displays the sky distribution of the lines of 
sight included in the statistical RC sample, demonstrat¬ 
ing the excellent coverage of the Galactic plane at low 
latitude. The fields are overlayed on the extinction map 


to a distance of = 5 kpc from Green et al. (2015) to il¬ 
lustrate the large amount of extinction affecting the disk 
region. While the distribution of fields includes a large 
number of fields in the region of the sky that contains the 
Galactic bulge, observations in these fields were limited 
to < 11. Because the RC has Mh ~ —1.49, this im¬ 
plies that the bulge is not actually included in our sample 
(these fields have D < 3kpc); the coverage is primarily 
of the disk. 

As discussed in detail in § 5 of Bovy et al. (2014), 
the overall-APOGEE and RC-specific selection cuts used 
to define the APOGEE-RC sample excludes stars with 
ages < 800 Myr and metallicities [Ee/H] < —0.9. The 
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Fig. 5.— Distribution of the 14,699 stars in the APOGEE-RC 
sample used in this paper in the plane defined by the iron abun¬ 
dance, [Fe/H], and the average abundance of a elements (see text). 
A linear binned density representation is used for the 68 % of 
the distribution that is contained within the shown contour. The 
dashed lines delineate the boundaries of the four broad subsamples 
that we study in § 4, which we denote with the given moniker. 

youngest part of the disk, as well as components more 
metal-poor than the bulk of the disk (e.g., Carollo et al. 
2010) are therefore not represented in the sample; how¬ 
ever, most of the disk mass satisfies these constraints. 
Because the RC consists of evolved giants stars that 
spend ^ 100 Myr of their lifetime in the RC phase, the 
population of RC stars is a biased tracer of the underly¬ 
ing stellar populations in the Milky Way. The extent of 
this bias is primarily dependent on the metallicity and 
the (unknown) age distribution, but the dominant effect 
is that the APOGEE-RC selection overrepresents stars 
in the 1 to 4 Gyr age range with respect to the underly¬ 
ing age distribution. Eor example, for a relatively con¬ 
stant star-formation history, the APOGEE-RG selection 
is expected to be dominated by stars with ages between 1 
and 4 Gyr, although older ages are represented at a lower 
level. 

Distances for stars in the APOGEE-RG catalog are 
determined by combining corrections to a single abso¬ 
lute magnitude as a function of ([J —i^s]o, [Ee/H]), 
based on stellar-evolution models (Bressan et al. 2012); 
these corrections are applied to each individual star. The 
overall distance scale is calibrated against an Hipparcos- 
based determination of the RG absolute Kg magnitude 
in the solar neighborhood (Laney et al. 2012). In this 
paper, we determine distances based on the i7-band lu¬ 
minosity rather than that in the i^g-band, because the 
APOGEE sample selection is performed in apparent H- 
band magnitude and our density-fitting formalism is sim¬ 
plest in this case. We follow the same procedure as dis¬ 
cussed in §§ 2 and 3 of Bovy et al. (2014) to compute the 
corrections to a single i^-band luminosity for the RG, as 
a function of {[J—Ks]o, [Ee/H])^, and to calibrate the dis¬ 
tance scale to the Hipparcos value of Mh = —1.49. The 

® We convert the metallicity [Fe/H] to metal mass fraction Z 


absolute magnitude Mh{[J—Ks]o, [Ee/H]) is displayed in 
Eigure 2. Eigure 3 displays the fractional difference be¬ 
tween the distances for stars in our sample determined 
from their Kg magnitude (the standard APOGEE-RG 
catalog product) and from their H magnitude. The frac¬ 
tional differences is < 3 % for all stars, with a median of 
2.3%. This is similar to the overall accuracy of the RG 
distance scale and so small that it does not impact the 
density fits below. 

2.2. Abundance measurements and main subsamples 

We make use of the abundance measurements provided 
in the DR12 release of APOGEE data (Holtzman et al. 
2015). Specifically, we use the iron content, [Ee/H], and 
the average abundance of a elements relative to iron. 
[Ee/H] is measured from 64 iron lines in the H’-band. As 
discussed by Holtzman et al. (2015), [Ee/H] is internally 
corrected for trends with Teff for stars in globular and 
open clusters observed by APOGEE. This internal cal¬ 
ibration only amounts to differences of ~ 0.04 dex over 
the ^ TOOK spanned by the RG. We determine the pre¬ 
cision in [Ee/H] by comparing measurements of [Ee/H] 
for different stars in the open clusters of Meszaros et al. 
(2013) in the Teff range of the RG and find that cr[Fe/H] = 
0.05 dex. 

We define the average a-enhancement as an average 
of the abundance of O, Mg, Si, S, and Ga. We do not 
include Ti, because of issues with its measurement in 
the H band (see Holtzman et al. 2015 for further discus¬ 
sion). Specifically, we average the abundances of [0/H], 
[Mg/H], [Si/H], [S/H], and [Ga/H] and subtract [Ee/H] 
to obtain the average ^-enhancement [a/Ee]. If no mea¬ 
surement was obtained for one of the five a elements, 
it is removed from the average. [0/H] is measured from 
the abundant molecular OH and GO features in the near- 
infrared that are, however, relatively weak for the warm 
RG stars. Mg, Si, S, and Ga are measured from neutral 
atomic lines for two (Ga) to 12 (Si) features. The abun¬ 
dances of a elements are similarly corrected for trends 
with Teff in clusters as [Ee/H] above; for the average 
[a/Ee] as defined here, this calibration only gives dif¬ 
ferences of 0.03dex over the ^ TOOK spanned by the 
RG. We determine the empirical precision in [a/Ee] us¬ 
ing the scatter in [a/Ee] for the calibration open clusters 
described above, and find that cr[Q,/Fe] ~ 0-02 dex, with a 
correlation of —0.4 with [Ee/H]. 

We apply a simple external calibration of [Ee/H] and 
[a/Fe] as follows. Using the APOGEE catalog abun¬ 
dances we determine the average [Ee/H] and [a/Ee] for 
giants observed by APOGEE in the temperature range 
of the RG (4500 K < Teff < 5200 K; see above) in 
M6T, an open cluster that provides an excellent chem¬ 
ical solar analog (e.g., Onehag et al. 2014). We find that 
[Ee/H]M 67 = 0.10 dex and [(a/Ee]M 67 = 0.03 dex; the 
large offset in [Ee/H]M 67 is at least partly due to an in¬ 
correct line list used in DR12 (see Shetrone et al. 2015). 
To calibrate APOGEE closer to the solar abundance 
scale using M6T, we apply constant offsets of —0.1 dex 
in [Ee/H] and —0.05dex in [a/Ee]. The reason for ap¬ 
plying an offset of —0.05 dex rather than —0.03 dex in 
[a/Ee] is that we define MAPs below using bins with 

assuming Zq = 0.017 and solar abundance ratios. 




7 






solar 

>: 


0 2 4 6 8 10 12 14 16 
R (kpc) 




low [Fe/H] •••. 

>; 


^0 2 4 6 8 10 12 14 16 

R (kpc) 


Fig. 6. — Observed spatial distribution of stars in the four broad subsamples defined in Figure 5, in Galactocentric radius, R, and distance 
from the mid-plane, Z. The Sun is located at (R, Z) = (8, 0.025) kpc. 


A[a/Fe] = 0.05 dex; a calibration offset of —0.05dex 
therefore does not change the binning of the data, which 
is the same when using the catalog abundances or the 
externally-calibrated abundances defined here. 

Figure 4 demonstrates the precision of our [tt/Fe] abun¬ 
dances. This figure considers RC stars with S/N > 200 
with —0.35 < [Fe/H] < —0.25. We interpolate the spec¬ 
tra for the 490 stars in our sample with these proper¬ 
ties onto a common Teff = 4800 K, log^ = 2.85 and 
[Fe/H] = -0.3, using quadratic interpolation. This re¬ 
moves the dominant effects of these stellar properties to 
reveal the effect of [a/Fe] on the spectra. We perform 
a PC A analysis of the residuals from a mean spectrum 
after the interpolation taking the measurement uncer¬ 
tainties into account (Bailey 2012), and only retain the 
components of each spectrum corresponding to the eight 
PC A eigenvectors with the largest eigenvalues. These 
eight eigenvectors explain all of the variance in the resid¬ 
uals above the measurement noise. Figure 4 displays 
median-stacked residuals of twelve random stars each in 
five bins in [a/Fe] with A[(a/Fe] = 0.05 dex ranging from 
[a/Fe] = 0.00 to [a/Fe] = 0.20. We perform this median 
stacking of a small number of spectra to reduce the noise 
for display purposes. It is clear from this figure that we 
measure [a/Fe] at very high precision from the spectral 
lines of O, Mg, Si, S, and Ca. The relative abundances 
of these elements with respect to each other appear to 
have little variation. 

Figure 5 displays the distribution of the statistical 
APOGEE-RC sampling in ([Ee/H], [a/Ee]). This dis¬ 
tribution contains the two main sequences—high- and 
low-[a/Ee]-seen in other investigations of this distribu¬ 
tion and discussed in detail for the APOGEE sample 
by Nidever et al. (2014) and Hayden et al. (2015). Also 
indicated in this figure are the four broad abundance- 
selected subsamples that we consider in § 4. The spatial 
distribution of stars in these four subsamples is displayed 
in Eigure 6. With due caution about the effect of selec¬ 
tion biases (which we correct for in the remainder of this 
paper) it is clear that the low-[a/Ee], low-[Ee/H] stars are 
primarily found in the outer disk, while the higher-[Ee/H] 
low-[a/Fe] stars are found closer to the center. High- 
[a/Fe] stars are found throughout the observed volume 
and are, in particular, more numerous at large distances 
from the mid-plane. 

In § 5, we consider MAPs, defined in the same manner 
as in Bovy et al. (2012c) as stars in bins with A [Fe/H] = 
0.1 dex and A[a/Fe] = 0.05 dex. From the discussion of 
the uncertainties above, it is evident that the contamina¬ 


tion between neighboring MAPs is slight, and that that 
between non-neighboring bins is essentially non-existent. 

3. DENSITY-FITTING METHODOLOGY 
3.1. Generalities 

To fit the spatial density profile for subsamples of RC 
stars, we follow the methodology of Bovy et al. (2012c) 
(see also Rix & Bovy 2013), who model the observed rate 
of stars in the joint parameter space of position, mag¬ 
nitude, color, and metallicity using a Poisson process. 
Best-fit parameters and their uncertainties for param¬ 
eterized spatial profiles are obtained by sampling this 
Poisson process’ likelihood of the observed data mul¬ 
tiplied with an uninformative prior. As discussed by 
Bovy et al. (2014) and above, for the RC stars we deter¬ 
mine the absolute magnitude Mh as a function of color 
(J — Ks)o and metallicity [Fe/H]; therefore, almost the 
same methodology for fitting the spatial profiles of RC 
stars applies here as was used by Bovy et al. (2012c) to 
analyze G dwarfs (whose absolute magnitude Mr was 
similarly determined from the color {g — r)o and [Fe/H]). 

In the present application, we therefore model the rate 
function X{O\0) that is a function of O = (/, 5, D, H, [J — 
i^s]o, [Fe/H]), and is parameterized by parameters 0; 
A(0|^) is given by 

X{O\0) = U,{X,Y,Z\0) X |J(X,T,Z;/,6,T))| 

X p{H, [J - Ks]o, [Fe/H] |X, T, Z) x S{1, 5, H ), 

( 1 ) 

where n^{'\0) is the spatial density in Galactocen¬ 
tric rectangular coordinates (A, T, Z) that we are ul¬ 
timately most interested in and that depends on 
parameters 6>, | J(A, T, Z;/, 6, T))| is the Jacobian of 
the transformation between (A, T, Z) and (/,5, U), 
p(i7, [J — As]o, [Fe/H] I A, T, Z) is the density of stars in 
magnitude-color-metallicity space, and 5'(/,6, iJ) is the 
survey selection function (the fraction of stars from the 
underlying population of potential targets observed by 
the survey). In APOGEE, the survey selection function 
is a function of position on the sky, is constant with (/, b) 
within an APOGEE field, and is a piecewise-constant 
function of apparent iJ-band magnitude, because tar¬ 
gets were sampled in magnitude bins. The APOGEE 
selection function is determined, tested, and discussed 
in detail in § 4 of Bovy et al. (2014). For the present 
work, we have updated the selection function to the full 
three-year data set presented in DR12. 

As in Bovy et al. (2012c), the rate has an additional 










































amplitude parameter. To remove this parameter from 
further cousideratiou, we margiualize the probability of 
the parameters of the rate fuuctiou over the amplitude 
of the rate. The margiualized likelihood cau be writteu 
as 


£(6») = \nv4Xi,Yi,Zi\e) -\nj d0A(0|6>) 


( 2 ) 


measuremeuts here. 

3.2. Stellar Number-density Models 

We fit a variety of models for the stellar uumber deusity 
of abuudauce-selected populatious. All of the models 
for which results are giveu iu this paper assume that 
the deusity is axisymmetric aud that the radial profile is 
separable from the vertical profile, such that 


Iu this expressiou, we have made use of the fact that 
the rate \{Oi\0) ouly depeuds ou 0 through aud 

therefore the Zi\0) is the ouly factor iu \{Pi\0) 

that depeuds ou 0] the other factors cau be dropped. 
The iutegral iu this equatiou is the effective volume of 
the survey that provides the uormalizatiou of the rate 
likelihood. It does uot depeud ou the iudividual data 
poiut, but is iustead a property of the whole survey for 
a giveu model specified by 0. 

Equatiou (2) is similar to the equivaleut likelihood iu 
Bovy et al. (2012c) (their equatiou [8]). The ouly sig- 
uificaut differeuce betweeu the two expressious is that 
the APOGEE selectiou fuuctiou depeuds ou the appar- 
eut maguitude H that is not corrected for extiuctiou, 
while the SEGUE selectiou was performed iu extiuctiou- 
corrected maguitude. Therefore, to calculate the uormal¬ 
izatiou iutegral iu equatiou (2) we require a model for the 
three-dimeusioual distributiou of extiuctiou Ah to cou- 
vert the model predictiou’s Hq to H. We discuss the 
methodology for efficieutly calculatiug the effective vol¬ 
ume for differeut types of surveys iu Bovy et al. (2016). 
Eor a peucil-beam survey like APOGEE, for which we 
cau assume that the deusity is coustaut over the area of 
each field, the effective volume cau be efficieutly com¬ 
puted as 


J d0A(0|6») = 

[ dDD‘^i^4[X,Y,Z]{D,\occition)\e) 

fields 


B(locatiou, D ), 


( 3 ) 


z/*(i?,0,Z) = E(i?)C(^), (5) 

where we defiue C{Z) such that J dZ C{Z) = 1. We have 
performed fits that add a coustaut deusity to capture auy 
outliers, but hud that the coutributiou from outliers is 
uegligible iu all cases; therefore we do uot iuclude them 
iu this descriptiou. The basic model for Y{R) that we 
cousider here is that of a brokeu expoueutial 


luE(i?) oc 


-^o) — -^peak 5 

-^0) R ^ -^peak • 


( 6 ) 


The relative uormalizatiou of these expoueutials is set 
to produce a coutiuuous Y{R) at i^peak- In MGMG ex- 
ploratious of the parameter coustraiuts, the iuverse scale 
leugths aud the logarithm of Rpeak are giveu flat priors. 
We also cousider models where Y{R) is a siugle expoueu¬ 
tial; such models esseutially have i^peak = 0. 

We have explored additioual fuuctioual forms for E(i?), 
such as a Gaussiau ceutered ou i?peak, but hud that the 
deusity is typically best described as a brokeu expoueu¬ 
tial. However, because the deusity drops quickly as oue 
moves away from i^peak, aud because the disk has a large 
rauge of i^peak (see below), determiuiug the exact form 
of the radial profile is difficult with the curreut data. 

We cousider four distiuct fuuctioual forms for C{Z). 
The first is that of a siugle expoueutial, with a scale 
height that is au expoueutial fuuctiou of radius R (i.e., 
a fiariug model) 


InC(Z) ^ exp [i? - i?o]) jZj . (7) 

To iuvestigate the fiariug profile further, we also cousider 
uou-expoueutial fiariug that is either liuear or iuverse- 
liuear, i.e. 


where 0(locatiou, D) is the effeetive seleetion funetion. 
Eor APOGEE, this is giveu by 

0(locatiou, D) = 5'(locatiou, k) 
k 

^(-^min,fe — Ho{D) < Ah{ 1, b, D) < H'max,/c ~ Ho{D)) 


where the sum is over the differeut maguitude bius, 
[^min,/c 7 ^max,/c]: fh^t stars are selected iu aloug each 
hue of sight. The uumerator U(i7niin,/c — < Ah{D) < 

Hmax,k—Ho) is the area of the APOGEE field iu questiou, 
with Ah betweeu the giveu bouudaries, aud the deuom- 
iuator Qf is the total area of the field. 5'(locatiou, k) is 
the APOGEE survey selectiou fuuctiou, i.e., the fractiou 
of poteutial targets with APOGEE spectroscopic obser- 
vatious iu each maguitude biu. This equatiou assumes 
that the RG is a staudard caudle with Mh = —1.49, al- 
lowiug us to compute Ho{D). Bovy et al. (2016) demou- 
strate that this assumptiou does uot affect the deusity 


InCiZ) c (1 ± [R - IZI , (8) 

where the sigus are such that is always uegative for 

outwardly-iucreasiug hz- In MGMG exploratious of the 
PDE, the iuverse scale height aud the iuverse fiariug 
scale leugth are giveu uuiform priors. A uou-fiariug 
vertical profile has = 0. The secoud form for ((Z) 

is a sum of two expoueutials 

^ 1^1) + exp (/iy 2 l^l) • 

(9) 

We also briefly cousider the geueralizatiou of the models 
iu equatious (7) aud (9), i.e., the two-expoueutial model 
where each scale height flares expoueutially as iu equa¬ 
tiou (7). As demoustrated below, the two-expoueutial 
form for the vertical profile does uot fit the data as well 
as the fiariug siugle-expoueutial form. Our maiu use of 
the secoud form is to iuvestigate whether auy MAPs show 
evideuce of a secoudary compoueut with a differeut scale 
height. 
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TABLE 1 

Results for broad abundance-selected samples 


Sample 

Density model 

Extinction map 

ft,in 
(kpc“^) 

, -1 
^ J?,out 
(kpc“^) 

^peak 

(kpc) 

hz 

(kpc) 

^flare ^2 

(kpc~ ^ or -) 

^Z,2 

(kpc) 

<1 

low [Fe/H] 

broken exp. 

w/ flare 

Marshall et al. (2006) 

0.27 ± 0.03 

0.36 ± 0.11 

10.8 ±0.1 

0.37 ± 0.02 

-0.09 ± 0.01 


0 

(low [«/Fe]) 



Green et al. (2015) 

0.28 

0.37 

10.8 

0.38 

-0.08 


13 




Sale et al. (2014) 

0.28 

0.37 

10.8 

0.38 

-0.08 


91 




Drimmel et al. (2003) 

0.27 

0.45 

10.9 

0.37 

-0.10 


649 




zero 

0.34 

0.59 

10.9 

0.39 

-0.10 


2162 


single exp. 


Marshall et al. (2006) 


0.06 ± 0.01 


0.45 ± 0.01 



723 


broken exp. 

w/ 2 hz 

Marshall et al. (2006) 

0.28 ± 0.03 

0.37 ± 0.02 

10.8 ±0.1 

0.47 ± 0.01 

< 0.01 

0.19 ± 97.44 

74 

solar 

broken exp. 

w/ flare 

Marshall et al. (2006) 

0.09 ± 0.04 

0.65 ± 0.02 

9.2 ±0.1 

0.28 ± 0.01 

-0.09 ± 0.01 


0 

(low [oi/Fe]) 



Green et al. (2015) 

0.15 

0.65 

9.2 

0.29 

-0.08 


77 




Sale et al. (2014) 

0.15 

0.65 

9.2 

0.29 

-0.08 


149 




Drimmel et al. (2003) 

0.07 

0.71 

9.4 

0.29 

-0.09 


286 




zero 

0.29 

0.78 

9.3 

0.32 

-0.07 


1443 


single exp. 


Marshall et al. (2006) 


0.33 ± 0.01 


0.31 ± 0.01 



906 


broken exp. 

w/ 2 hz 

Marshall et al. (2006) 

0.11 ± 0.03 

0.66 ± 0.02 

9.2 ±0.1 

0.33 ± 0.01 

< 0.01 

0.11 ± 29.94 

90 

high [Fe/H] 

broken exp. 

w/ flare 

Marshall et al. (2006) 

0.28 ± 0.15 

0.81 ± 0.03 

6.8 ± 0.2 

0.27 ± 0.01 

-0.14 ± 0.02 


0 

(low [oi/Fe]) 



Green et al. (2015) 

0.62 

0.79 

6.5 

0.28 

-0.12 


98 




Sale et al. (2014) 

0.56 

0.79 

6.6 

0.28 

-0.12 


148 




Drimmel et al. (2003) 

0.20 

0.81 

6.8 

0.28 

-0.14 


84 




zero 

1.07 

0.78 

6.5 

0.30 

-0.10 


832 


single exp. 


Marshall et al. (2006) 


0.60 ± 0.01 


0.28 ± 0.01 



430 


broken exp. 

w/ 2 hz 

Marshall et al. (2006) 

0.43 ± 0.32 

0.80 ± 0.08 

6.6 ± 1.7 

0.28 ± 0.01 

< 0.02 

0.09 ± 1.66 

95 

high [a/Fe] 

broken exp. 

w/ flare 

Marshall et al. (2006) 

0.93 ± 0.57 

0.43 ± 0.03 

< 4.4 

0.95 ± 0.05 

0.01 ± 0.03 


0 




Green et al. (2015) 

1.16 

0.43 

1.0 

0.97 

0.02 


6 




Sale et al. (2014) 

0.95 

0.43 

1.8 

0.97 

0.02 


10 




Drimmel et al. (2003) 

0.88 

0.45 

4.1 

0.96 

0.00 


14 




zero 

0.00 

0.40 

1.1 

1.09 

0.04 


170 


single exp. 


Marshall et al. (2006) 


0.43 ± 0.02 


0.95 ± 0.05 



0 


broken exp. 

w/ 2 hz 

Marshall et al. (2006) 

0.90 ± 0.57 

0.43 ± 0.02 

2.2 ± 1.2 

0.95 ± 0.06 

< 0.02 

0.15 ± 195.28 

0 


Note. — Lower limits are at 95% posterior confidence. The seventh column is ^^e flaring model and the amplitude /32 for the model with two vertical 

scale heights. The parameters /32 and h,^ 2 are marginalized under the constraint that h,^ 2 ^0 9o different from (to avoid the massive degeneracy when 

and 2 3-*^® allowed to be the same). The model consisting of a broken radial exponential with a flaring exponential scale height provides the best fit in all 
cases, although the high-[Q:/Fe] sample is always fit as a single radial exponential. 


We optimize the likelihood in equation (2) for each 
density model and each data set below using a downhill 
simplex algorithm. We then use this optimal solution to 
initiate an MCMC sampling of the posterior PDF, ob¬ 
tained using an affine-invariant ensemble MCMC sampler 
(Goodman & Weare 2010; Foreman-Mackey et ah 2013). 
Reported parameter estimates are based on the median 
and standard deviation of one-dimensional projections of 
the MCMC chain. 

3.3. Tests on mock data 

We have performed a suite of mock-data tests to vali¬ 
date our code (i.e., checking that we recover the correct 
input density profiles when fitting a model that includes 
the input assumptions) and to determine the impact of 
the uncertainty in the three-dimensional extinction map. 
We generate mock data with profiles similar to a few of 
the best-fit profiles for the broad abundance-selected sub¬ 
samples below. In particular, we produce mock data with 
a fiat radial profile, a single-exponential profile with a 
scale length of 3 kpc, or broken-exponential profiles with 
peak radii of 8 and 11 kpc; all mock data have a con¬ 
stant thickness with a scale height of 300 pc, except that 
we also generate mock data with the exponential flaring 
profile of equation (7) with a scale length of 10 kpc for 
the broken-exponential profiles. All mock data assume 
extinction according to the Green et al. (2015) map and 
properly take into account the variation of the density 
profile within each APOGEE field, therefore testing that 
this can be ignored in the calculation of the effective 
volume (see above). Each mock-data sample has 20,000 
stars. 

When fitting the mock data with profiles that include 
the input density profile, we find that our code recov¬ 
ers the correct profile to within the uncertainties that we 
find for the real data below. In particular, we confirm 
that radial profiles consisting of a single exponential are 


recovered as such, even when fit with more the general 
broken-exponential profile. Eurt her more, we recover that 
no constant-scale-height mock data set displays any flar¬ 
ing and that the flaring scale length is correctly recovered 
in the mocks with flaring. Eor the latter, we also find that 
the radial profile is almost perfectly recovered even when 
they are fit with a constant scale height. This confirms 
the expectation that the radial and vertical profiles are 
measured almost independently of each other due to the 
many lines of sight at 6 = 0° in the APOGEE footprint 
(see Eigure 1). 

To determine the impact of the uncertainty in the 
interstellar extinction, especially in the inner MW, we 
fit the mock data assuming extinction according to 
Marshall et al. (2006) near the mid-plane, rather than 
the Green et al. (2015) extinction assumed in generating 
the mock data; the latter is typically smaller than the 
former. We find only small changes in the best-fit pa¬ 
rameters when using the incorrect extinction map; these 
differences are smaller or similar to the uncertainties in 
the best-fit parameters, and do not lead to qualitatively- 
different inferred density profiles. Even employing the 
map of Drimmel et al. (2003) in the fit only changes the 
best-fit parameters by insignificant amounts, similar to 
what we find for the real data below. Therefore, we con¬ 
clude that the uncertainty in the three-dimensional ex¬ 
tinction does not limit our understanding of the spatial 
density profiles of stellar subsamples in APOGEE. 

4. THE SPATIAL STRUCTURE OF BROAD 
ABUNDANCE-SELECTED SUBSAMPLES 

We now discuss the results from fitting the spatial den¬ 
sity of the broad abundance-selected subsamples (Eigures 
5 and 6), using the functional forms described in § 3.2 
above. We start by fitting broad sub-samples first, be¬ 
cause of the (relatively) large sample size that allows us 
to precisely determine the shape of their density profiles. 
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Fig. 7.— Comparison between the observed (filled histogram) distribution of distance moduli, /r, and that predicted by the best-fit models 
for the three low-[a/Fe] subsamples of Figure 5 (rows) and different regions of the sky (columns). The red curves demonstrate the fits 
using a combination of the Marshall et al. (2006) and Green et al. (2015) 3D extinction maps for the different density profiles displayed 
in Table 1. The remaining curves display the fits assuming different 3D extinction maps. The predicted number of stars in each spatial 
region is given for the fiducial model (broken exp. w/ fiare with Marshall et al. 2006 extinction). Different 3D extinction maps give by 
and large similar fits, except for in the inner Galaxy (left panels), where the Marshall et al. (2006) map performs best. A model with zero 
extinction (black curves) provides poor fits for all low Galactic latitude locations. The sharp features in the zero-extinction model refiect 
the discontinuous nature of the APOGEE selection function as a function of H; these are smoothed by the extinction for models with 
extinction (see Figure 6 in Bovy et al. 2016). The density in all cases is best fit as a radially broken exponential with a fiaring vertical scale 
height; a single radial exponential fails to explain the radial profile over all radii, demonstrated by the poor fit in the outer regions of the 
disk (third panels). 


In § 5 we then refine our results by fitting the preferred 
models from this section to the MAPs, which have much 
smaller sample sizes. The well-populated broad subsam¬ 
ples also make it possible to demonstrate the goodness of 
the fits of different spatial profiles by directly overlaying 
the best-fit models on the observed distribution of stars. 

4.1. The surface-density profile 

Results from fitting the density profiles given in § 3.2 
to the broad abundance-selected subsamples are given in 
Table 1. Comparisons between different model fits dis¬ 
cussed below and the observed distance distribution are 
displayed in Figure 7 for the three low-[(a/Fe] subsam¬ 
ples, and in Figure 8 for the high-[(a/Fe] subsample. We 
fit the radial and vertical profiles simultaneously, but fo¬ 


cus on discussing the resulting radial profiles S(R) in this 
subsection. 

It is clear from Table 1 and Figure 7 that a single¬ 
exponential radial profile S(R) (dotted line in Figure 7) 
does not provide a good fit to the data for the low-[(a/Fe] 
subsamples. A broken exponential of the type as in 
equation (6) provides a much better fit. Discrepancies 
between the observed and predicted distribution of dis¬ 
tances in Figure 7 are small for this model. In addition 
to the broken-exponential model, we have also performed 
fits with a Gaussian radial profile. These only gave sim¬ 
ilar or worse fits to the data, but with of only tens. 
The exact shape of the radial profile can therefore not be 
determined at high confidence, but it is clear that S(R) 
for the low [a/Fe] subsamples increases up to a peak 
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Fig. 8. — Same as Figure 7, but for the high-[Q;/Fe] sample. All different density profiles provide equally good fits, demonstrating that 
the density is very close to a single radial and vertical exponential for this population. 



Fig. 9.— Radial surface profile 51(77) of the four broad abundance-selected subsamples indicated in Figure 5. The gray region gives the 
95 % uncertainty range. All profiles are relative to the density at 77 = 8kpc; an arbitrary offset in the vertical direction has been applied to 
separate the four profiles. The three low-[Q;/Fe] subsamples are best represented as a broken exponential, while the high-[Q;/Fe] subsample 
consists of a single exponential distribution over the full radial range that is observed. 


radius i^peak and declines beyond that. S(i?) is highly 
inconsistent with being a single exponential. 

For the high-[(a/Fe] subsample the single-exponential 
model provides a good fit, while fits with a broken- 
exponential or other radial profile all end up as close 
to a single exponential as possible. For example, the 
broken-exponential fit places i^peak outside of the ob¬ 
served volume (i^peak < 4.4 kpc; Table 1) such that this 
fit is equivalent to a single exponential. The Gaussian ra¬ 
dial profile does the same, and adjusts the width param¬ 
eter such that the profile closely approximates a single 
exponential. We can therefore be confident that S(i?) 
for the high-[a/Fe] subsample is very close to a single 
exponential. 

The fits in Table 1 and Figures 7 and 8 are performed 
for several different extinction maps. The standard ex¬ 
tinction map used is that of Green et al. (2015). Used on 
its own it is labeled as “Green et al. (2015)”. When we 
replace the part of it at —100° < I < 100° and |6| < 10° 
with the map of Marshall et al. (2006), we label this as 
“Marshall et al. (2006)”; Likewise, when we replace the 
part of it that overlaps with the map of Sale et al. (2014), 
we refer to that model as “Sale et al. (2014)”. We also 
test the performance of two other extinction maps: that 
of Drimmel et al. (2003), which is defined over the entire 


sky, and a model without any extinction (labeled “zero”). 

In all cases the combination of the Marshall et al. 
(2006) and the Green et al. (2015) maps provides the 
best fit. The map of Sale et al. (2014) performs slightly 
worse than that of Green et al. (2015) where Sale et al. 
(2014) overlaps the latter. The map of Drimmel et al. 
(2003), which consists of a simple model for the three- 
dimensional distribution of dust and has lower angular 
resolution than the other maps, gives very poor fits. The 
model with zero extinction clearly provides a bad fit to 
the data, both from the Ax^ in Table 1, and directly 
from the comparison between the model and the data in 
Figure 7, especially at I < 70°. All of the different extinc¬ 
tion maps, however, give very similar best-fit parameters 
for the basic model for the spatial density. 

The radial profile and its uncertainties for the stan¬ 
dard model of a broken exponential for the four broad 
subsamples is displayed in Figure 9. It is clear that the 
uncertainty in this (parametric) model is small, and that 
the peak radius i^peak for the low-[(a/Fe] subsamples is 
larger for decreasing [Fe/H]. The shape of the radial pro¬ 
file around i^peak is quite similar for all three low-[(a/Fe] 
subsamples, with a shallow rise at i? < i^peak and a steep 
decline at > i^peak- The high-[a/Fe] subsample could 
be thought of as having i?peak ^ 4 kpc, and therefore be 
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the continuation of the trend of the low-[(a/Fe] subsam- 
ples, but with the current radial coverage we cannot test 
that scenario. 


4.2. The vertical profile 

Having determined that a broken-exponential T^{R) 
provides the best fit, we fit three different vertical pro¬ 
files to the APOGEE-RC data. The simplest model is 
that of a single-exponential model with a radially con¬ 
stant scale height hz- By fitting more complex models, 
we find that this model is strongly ruled out for the low- 
[a/Ee] subsamples. They are instead better fit with a 
model where hz is a function of R, and we employ the 
flaring model of equation (7) (the flaring profile is ex¬ 
plored in more detail for the MAPs below using the al¬ 
ternative flaring models). An alternative model is that 
each subsample consists of the sum of two exponentials 
(the model of equation [9]); this model does not fit as 
well (see Table 1). We have also fit a generalized model 
where the vertical profile consists of the sum of two ex¬ 
ponentials that flare exponentially with the same scale 
length. In all cases, the best fit for this general model 
reverts to that of the single-exponential, flaring model. 
All of the low-[a/Ee] subsamples are consistent with a 
common flaring scale length of ~ —O.lkpc”^. We 

refine this measurement in § 5.2 below. 

Like for the surface-density profile, the high-[(a/Ee] 
subsample is consistent with the simplest model, in this 
case a single vertical exponential with a constant hz{R)- 
That is, we see with high confidence that the high-[Qf/Ee] 
subsample does not display the same kind of flaring as the 
low-[(a/Ee] subsamples, but hz{R) is instead constant. 
We refine the quantitative constraint in § 5.2 below. 

5. THE SPATIAL STRUCTURE OF MAPS 

In this section we repeat the density fits from the previ¬ 
ous section, but we perform them on abundance-selected 
subsamples that are narrower in [Ee/H] and [a/Ee]. That 

is, we use MAPs, defined here as abundance bins with 
widths of A[Ee/H] = 0.1 dex and A[a/Fe] = 0.05 dex. 
We do not make use of other abundances for defining 
MAPs, but stress that our empirical description does not 
assume chemical homogeneity beyond ([Ee/H], [a/Ee]). 
As discussed at the end of § 2.2, these widths are about 
twice as large as the uncertainties in these quantities 
and the contamination between MAPs is therefore small. 
Because of the small number of stars in the statistical 
APOGEE-RC at high [a/Ee], we perform fits for MAPs 
with at least 15 stars; the measurements for MAPs with 
so few stars are noisy, but informative enough to help 
establish trends. We again discuss the results for the 
surface-density and vertical profiles separately, but both 
were measured simultaneously for all MAPs. 

5.1. The surface-density profile 

Inspired by the fits to the broad abundance-selected 
subsamples in § 4.1, we fit a broken-exponential E(R) to 
each well-populated MAP. We constrain these broken- 
exponential models to have an inner profile that is in¬ 
creasing with R and an outer profile that is decreas¬ 
ing, although the vast majority of MAPs have well con¬ 
strained profiles without this prior constraint that satisfy 

it. Eor MAPs that are best fit as a single exponential. 



-0.6 -0.4 -0.2 0.0 0.2 


|Fe/H] 

Fig. 10. — Peak radius of the radial profiles of MAPs. This figure 
displays the Galactocentric radius at which the surface density of 
each MAP peaks when fit with a broken-exponential radial profile 
that is constrained to be increasing within Rpeak- The locus of 
the high- and low-[a/Fe] sequences in Figure 5 are indicated with 
gray curves. MAPs along the high-[a/Fe] sequence do not display 
a peak in their radial profile within the observed R range; they are 
best represented as a single radial exponential and are indicated 
as “Rpeak < 5kpc”. MAPs along the low-[a/Fe] sequence show a 
striking increase in Rpeak with decreasing [Fe/H]. 


this constraint forces Rpeak to lie at small R; without 
this constraint, the fit would have a degeneracy between 
very small Rpeak and very large Rpeak (as long as Rpeak 
is outside of the observed volume). We always use the 
combined Marshall et al. (2006) and Green et al. (2015) 
extinction map, which provided the best fit to the broad 
subsamples above. 

We display the dependence on ([Ee/H], [<a/Ee]) of the 
peak of E(R) in Eigure 10. We determine Rpeak typically 
to 0.3 kpc, while the range of Rpeak covers about 8kpc. 
Thus, the smooth trends seen in Eigure 10 are deter¬ 
mined at high significance. These more detailed results 
confirm the behavior found for the broad subsamples in 
§4.1. Low-[(a/Ee] MAPs have an increasing Rpeak with 
decreasing [Ee/H], ranging from Rpeak ~ 5kpc at the 
highest [Ee/H], to Rpeak ~ 13kpc at the lowest [Ee/H]. 
This trend has a weak dependence on [a/Ee]. We have 
indicated the locus where the low-[a/Ee] sequence is well 
populated (that is, where Rpeak is best determined) and 
along this sequence the correlation between [Ee/H] and 
Rpeak is incredibly tight. 

The behavior of individual high-[(a/Ee] MAPs also con¬ 
firms that high-[(a/Ee] MAPs do not display a break in 
their surface-density profiles, but are instead consistent 
with a single exponential. Eor all MAPs along the high- 
[a/Ee] sequence, Rpeak is constrained to lie outside of the 
observed volume. 

The radial dependence of E(R) is displayed in Eig¬ 
ure 11. This figure only shows MAPs along the well- 
populated high- and low-[(a/Ee] sequences for clarity, but 
the behavior of other MAPs is similar, albeit noisier. It is 
clear that the radial profiles for all but the lowest [Ee/H] 
MAP are well constrained to have the broken-exponential 
form with an almost universal shape around Rpeak- The 
inner profile is typically shallower than the outer profile, 
except for lower-[Ee/H] MAPs. However, these MAPs 
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Fig. 11.— Radial surface profile of MAPs. For display purposes, MAPs along the well-populated low- and high-[Q;/Fe] sequences are 
shown in the bottom and top panel, respectively, but the trends are the same for all MAPs. For the low-[Q;/Fe] sequence these are the 
MAPs with [a/Fe] = +0.05 up to [Fe/H] = —0.4 and [a/Fe] = 0.0 at higher [Fe/Hj. For the high-[a/Fe] sequence these are the MAPs with 
[a/Fe] = +0.20 for [Fe/H] = (-0.5,-0.4), [a/Fe] = +0.15 for [Fe/H] = (-0.3,-0.2), and [a/Fe] = +0.10 for [Fe/H] = -0.1. The colored 
bands give the 95 % uncertainty region. The radial profiles of high-[a/Fe] MAPs are given by single exponentials with a common scale length 
of 2.2 ± 0.2 kpc. The metal-poor low-[a/Fe] MAPs peak in the outer disk (Rpeak ^ lOkpc) and are spread over a wide range of radii, with 
a relatively shallow outer scale length of about 3kpc. The metal-rich low-[a/Fe] MAPs are very centrally-concentrated (Rpeak 8kpc), 
with outer scale lengths of only ps 1.25 kpc. The inner, rising scale length is universally 3 kpc, in all MAPs where it can be constrained. 


are only sparsely populated at the large distance from 
the Galactic center where their outer profiles are con¬ 
strained, as is also evident from the uncertainties. The 
top panel displays S(i?) for the high-[a/Fe] MAPs. In 
addition to all being consistent with a single exponential, 
they are all consistent with the same single exponential. 


with a scale length of Hr = 2.2 ± 0.2 kpc. 

5.2. The vertieal profile 

Results from fitting the MAPs with the standard model 
for the vertical density—a single exponential with an 
exponentially-increasing scale height—are displayed in 
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Fig. 12. — Flaring of individual MAP disk components. The 
PDFs for the inverse flaring scale length for individual MAPs 

for well-populated high-[Q:/Fe] (top panel) and low-[a/Fe] (bot¬ 
tom panel) MAPs are displayed (histogram), together with smooth 
fits to these PDFs with a sum of two Gaussians (colored curves). 
The combined PDF obtained by multiplying the individual smooth 
PDFs is shown in black; its mean and standard deviation are given 
in the top right. The displayed MAPs are those of Figure 11, with 
colors ranging from blue to red for low- to high-[Fe/H] MAPs. The 
dashed vertical line indicates the limit of flaring disks (to its left); 
the high-[a/Fe] MAPs display no flaring to high precision, while the 
low-[a/Fe] MAPs are consistent with a single flaring scale length 
of i^flare ~ ± 0.7 kpc. 

Figures 12, 13, 14, and 15. The vertical profile of individ¬ 
ual MAPs is difficult to determine from the APOGEE- 
RC data, because the sample is dominated by low- 


latitude fields that give little leverage for measuring the 
vertical density drop off. When allowing the inverse scale 
length of the flare to be free, there is a large 

degeneracy between scale height at the 

solar circle Eigure 12 displays PDEs for for 

MAPs along the well-populated high- and low-[a/Fe] se¬ 
quences. Eor these MAPs, relatively well con¬ 

strained. We see from these PDFs that the high-[(a/Fe] 
MAPs are all consistent with having a constant hz{R), 
while low-[(a/Fe] MAPs show strong evidence for a flaring 
hz{R)- The constraints on both [a/Fe] groups 

are tight: ~ ^ 0.02 kpc“^ for the high-[a/Fe] 

MAPs and R^l^^ = —0.12±0.01 kpc~^ for the low-[(a/Fe] 
MAPs. The latter corresponds to a flaring scale length of 
8.5 ± 0.7 kpc. We illustrate the flaring (and non-flaring) 
of the MAPs in Figure 13. To further explore what the 
radial dependence of hz is, we have also fit each MAP 
with the linear and inverse-linear flaring profiles of equa¬ 
tion (8). Figure 13 contains these results: the median 
of the MCMC samples from these fits are very close to 
that using the exponential model. The gray band for 
each MAP displays the range in hz spanned by the dif¬ 
ferent models. This gray band is only visible for a few 
low-[(a/Fe] MAPs and mainly near the edges of the ra¬ 
dial range, indicating that the low-[a/Fe] MAPs have an 
exponential flaring profile to within their uncertainties. 
For high-[(a/Fe] MAPs the results from the alternative 
flaring profiles are almost indistinguishable from those 
with the exponential model. 

To determine the scale heights, of all of the MAPs, 
we then repeat the density fits while fixing Specif¬ 
ically, we set = 0 for all of the high-[(a/Fe] MAPs, 

defined here as those with Rpeak < 5 kpc (see Fig¬ 
ure 10), and Rflare ~ —0.1kpc“^ otherwise. The result¬ 
ing hz([Fe/H], [<a/Fe]) is displayed in Figure 14. In addi¬ 
tion to the requirement that each MAP needs to have at 
least 15 stars, we have further removed MAPs for which 
the uncertainty in hz is larger than 20 %; this removes a 
few MAPs at low [Fe/H] (both at high and low [a/Fe]). 
It is clear from this figure that the dependence of hz 
on ([Fe/H], [a/Fe]) is very smooth, and similar to that 
found by Bovy et al. (2012c). In particular, intermedi¬ 
ate components with hz ~ 500 pc are prevalent at the 
low-[Fe/H] end of the low-[(a/Fe] sequence and at the 
high-[Fe/H] end of the high-[a/Fe] sequence. 

In § 4.2 we determined that a vertical profile consisting 
of the sum of two exponentials does not provide a good fit 
to the broad abundance subsamples. We repeat this ex¬ 
ercise here, fitting a model with two vertical exponentials 
with a constant scale height. In Figure 15 we compare 
the scale height of the dominant component with that of 
the secondary component. We constrain the secondary 
component to contribute at least 15 % of the local col¬ 
umn density, to remove the large degeneracies at small 
amplitudes of the secondary component. We find that 
whenever the secondary component makes a substantial 
contribution to the column density, its scale height is the 
same as that of the dominant component. That is, the 
vertical profile is consistent with being a single exponen¬ 
tial. 


6. DISCUSSION 
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Fig. 13.— Vertical profile of MAPs. This figure displays the radial dependence of the scale height of MAPs along the well-populated 
high-[Q:/Fe] (top panel) and low-[a/Fe] sequence (bottom panel] those of Figure 11, except for the lowest-[Fe/H] MAP). The colored bands 
give the 95% uncertainty region around the median displayed as a solid line for the exponential flaring model. The gray band (when visible) 
shows the range in median profile spanned by all three flaring models (exponential, linear, and inverse-linear; see equations [7] and [8]). 
A [Fe/H]-dependent offset has been applied to the y axis to separate the different MAPs; the dashed horizontal line gives the position of 
hz = 300 pc for each MAP. The high-[Q;/Fe] MAPs do not flare, while all low-[a/Fe] MAPs are consistent with an exponential flaring profile 
with a scale length of T^flare ~ 0.7kpc (see Figure 12 above). The scale height increases smoothly from high- to low-[Fe/H] MAPs. 

This paper is the first to dissect in detail—and over a extent to which these results are compatible with earlier 
wide range of Galactocentric radii—the radial structure analyses, before proceeding to interpret them in a galaxy 
of the MW’s stellar disk in terms of MAPs, abundance- formation context, 
selected stellar sub-populations. As some of the results of 
this analysis may appear unorthodox, we first show the 
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Fig. 14. — Vertical scale height of MAPs. The scale heights, hz^ 
for fits assuming a constant scale height with Galactocentric ra¬ 
dius for MAPs best fit as a single radial exponential (those labeled 
“i^peak < 5 kpc” in Figure 10) and for fits with a flaring vertical 
profile with a flare scale length of 10 kpc for all other MAPs. A 
few MAPs along the upper and left edges for which the scale height 
cannot be precisely determined from the present data are omitted. 
The locus of the high- and low-[a/Fe] sequences are indicated as in 
Figure 10. MAPs displays a smooth increase in hz as a function 
of declining [Fe/H] and increasing [a/Fe], even for the conservative 
flaring model chosen here. 



Fig. 15. — Scale height of the dominant component versus that 
of the secondary component in two-vertical-exponential fits to the 
MAP spatial densities. The scale heights displayed here are the me¬ 
dian of samples from the posterior PDF for which the secondary 
component provides significantly to the density (defined as con¬ 
tributing more than 15% of the local column density), to avoid 
the massive degeneracy in the scale height when the amplitude is 
allowed to be close to zero. The two scale heights are the same for 
all MAPs, demonstrating that they are well represented as a single 
vertical exponential. MAPs with very low scale heights, however, 
show some evidence of a second, even lower scale height. 

6.1. Comparison to Bovy et al (2012c) 

For the high-[a/Fe] MAPs our results on the radial 
profile are in perfect agreement with those of Bovy et al. 


(2012c). They find for a single-exponential fit to S(i?) for 
high-[(a/Fe] MAPs that Hr = 2.01 ± 0.05 kpc, consistent 
with our measurement here of Hr = 2.2 ± 0.2 kpc. Their 
result is more precise because they have ^ 14, 000 high- 
[a/Fe] stars compared to only 500 here. However, the 
RC distance scale is more accurately known than that 
of the G dwarfs used by Bovy et al. (2012c), so larger 
samples of RC stars should eventually lead to a more 
accurate and precise measurement of Hr. The improved 
radial coverage of the APOGEE-RC sample over the G- 
dwarf sample also allows us to ascertain that E(R) is 
indeed a single exponential; this was an assumption in 
Bovy et al. (2012c). 

The present analysis has shown that the radial surface- 
density profile of low-[(a/Ee] MAPs is not a single expo¬ 
nential, but is much better fit as a broken-exponential 
profile, rising to a peak radius Rpeak, before falling off. 
This makes it difficult to meaningfully compare our mea¬ 
surements of E(R) to measurements in the literature, 
which typically consist of single-exponential fits to a mix 
of MAPs. 

Bovy et al. (2012c) were the first to dissect the MW 
disk into narrow abundance bins and while they also fit 
single-exponential models to E(R), their results can be 
more easily compared to the current measurements. Eor 
comparison, we have carried out single-exponential fits to 
the APOGEE data E(R) (see Table 1), but the different 
radial coverage makes this a qualitative exercise rather 
than a quantitative one. On this level, our results are 
consistent with those of Bovy et al. (2012c) for the low- 
[a/Fe] MAPs. When fit with a single-exponential E(R), 
low-[Ee/H] MAPs have a flat profile, solar-[Ee/H] MAPs 
have a scale length of Hr = 3.0±0.1 kpc, and high-[Ee/H] 
MAPs have Hr = 1.67±0.03kpc, similar to the results in 
Eigure 5 of Bovy et al. (2012c). It is clear, however, that 
we significantly refine the results of Bovy et al. (2012c) 
for the low-[a/Ee] MAPs here by determining the shape 
of E(R). 

We find that the vertical profile of the low-[a/Ee] 
MAPs consists of a single exponential, but with a scale 
height that is flaring outward with an approximately ex¬ 
ponential profile. Bovy et al. (2012c) assumed a con¬ 
stant hz{R)^ because (unpublished) investigations using 
the SEGUE G-dwarf data set demonstrated that all but 
the most extreme flaring models were consistent with 
the data, and hz{R) = constant was the simplest possi¬ 
ble assumption. The inability to determine the flaring of 
low-[a/Ee] MAPs by Bovy et al. (2012c) was due to the 
limited radial coverage (spanning only about ±2 kpc) 
and the lack of low latitude lines of sight in SEGUE. We 
have not attempted to refit the SEGUE G-dwarf data 
using the best-fit flaring model from this paper, but the 
slow exponential flaring of = —0.12 ± 0.01 kpc~^ 

is such that it is most likely consistent with the SEGUE 
G-dwarf data. 

The measurements of the scale heights in this pa¬ 
per are much more uncertain than those of Bovy et al. 
(2012c), because most of the APOGEE lines of sight are 
concentrated near the plane. Because we do not cur¬ 
rently possess the relative calibration of ([Ee/H], [a/Ee]) 
in APOGEE and SEGUE we cannot directly compare 
measurements of hz in SEGUE and in APOGEE. How¬ 
ever, it is clear that the main trends are the same: 
























17 


(a) /iz([Fe/H], [a/Fe]) spans the range from 200pc to 
1000 pc, with a smooth transition between these two ex¬ 
tremes, (b) the thinnest components have low [a/Fe] and 
high [Fe/H], (c) the thickest components have high [a/Fe] 
and low [Fe/H], and (d) intermediate {hz ~ 500pc) com¬ 
ponents lie at the high-[Fe/H] end of the high-[(a/Fe] 
sequence and the low-[Fe/H] end of the low-[(a/Fe] se¬ 
quence. Point (a) here is particularly remarkable, be¬ 
cause we measure hz much closer to the plane and for 
isothermal MAPs (Bovy et al. 2012b) in dynamical equi¬ 
librium, we expect hz to become larger near the plane. 
We also confirm that each MAP consists of a single ver¬ 
tical exponential component. 

6.2. Comparison to star-count measurements 

Any quantitative comparison with previous measure¬ 
ments of S(i?) is problematic, because these measure¬ 
ments typically fit single-exponential S(i?), while we 
find that the radial profile of MAPs is better fit as a 
broken exponential with a peak radius. We are there¬ 
fore forced to qualitatively interpret older measurements 
based on their radial and vertical coverage. The most 
reliable measurements of the radial structure are based 
on IR data. Among these, GLIMPSE measurements 
based on |5| < 1° at 10° < I < 65° give hR = 
3.9 ± 0.6kpc (Benjamin et al. 2005), while COBE near- 
IR and far-IR disk data at |6| < 30° lead to hn = 2.3 kpc 
(Drimmel & Spergel 2001), and near-IR 2MASS data in 
the outer disk are best fit with hR = 2.2 kpc (Reyle et al. 
2009). We can explain these differences qualitatively 
as follows: the GLIMPSE star counts at \b\ < 1° are 
dominated by the low-[(a/Ee] MAPs, which in the inner 
MW are approximately a combination of the “solar” and 
“high-[Ee/H]” populations in Eigure 9 and therefore have 
quite fiat E(R). The COBE measurements that extend 
to higher \b\ contain a much higher contribution from 
the short Hr high-[a/Ee] MAPs, and therefore lead to 
hR closer to 2 kpc. The outer-disk 2MASS sample from 
Reyle et al. (2009) is dominated by the steep outer scale 
length of low-[(a/Fe] “solar” populations {hR ^ 1.5 kpc), 
and the somewhat steep outer profile of “low-[Ee/H]” 
populations {hR ^ 2.7kpc; see Table 1); that an overall 
profile with hR = 2.2 kpc would result appears somewhat 
likely. While these qualitative comparisons are interest¬ 
ing, they mostly point toward the necessity to perform 
proper comparisons between the radial profile obtained 
from different types of tracers and from different parts of 
the MW disk. 

Star counts in the outer disk have found evidence of a 
break in the stellar surface density (i.e., an abrupt steep¬ 
ening of the already decreasing surface density), with 
the majority of studies converging on a break around 
R ^ 13.5 kpc (e.g., Reyle et al. 2009; Sale et al. 2010; 
Minniti et al. 2011; R. Benjamin et ah, 2016, in prepa¬ 
ration), with a steep density decline beyond the break 
{hR = 1.2 ± 0.3kpc; Sale et al. 2010). We have found 
that each MAP has a peak radius, similar in kind to the 
break radius found in the outer disk, with outer scale 
lengths hR^out ^ 2 kpc (see Figure 11). The break radius 
seen in star counts using data and methods that do not 
discern stars by their abundances is therefore most likely 
nothing more than the outermost of the series of break 
(or peak) radii displayed by the MAPs; the series in Fig¬ 
ure 10 ends around 13 kpc. More metal-poor MAPs fill-in 


the density beyond the peak of more metal-rich MAPs 
until the most metal-poor MAPs are reached, and the to¬ 
tal surface-density starts to decline with the steep outer 
scale length of the most metal-poor low-[(a/Fe] MAPs. 
Because we simultaneously fit for E(R) and the flaring 
of the disk, we can be sure that the steep decline in mid¬ 
plane density is truly because of a declining E(R) and 
not just a consequence of the flaring of the disk. 

6.3. Implications for disk formation and evolution 

Our new results on the stellar population dependent 
structure of the MW disk provide stringent and qualita¬ 
tively new constraints on models of the formation and 
evolution of galactic disks. We discuss qualitatively the 
implications of our results for some of the main evolu¬ 
tionary mechanisms, but more detailed comparisons to 
models are beyond the scope of this paper. 

As displayed in Figures 10 and 11, we And that when 
the MW disk is dissected into MAPs, it consists of a set 
of donut-like rings, with a peak radius that is a declin¬ 
ing function of [Fe/H] for low-[(a/Fe] MAPs. high-[(a/Fe] 
MAPs have peak radii constrained to be < 5 kpc, and are 
therefore consistent with a more traditional disk struc¬ 
ture. Whether or not they actually display at peak at 
0 < R/kpc < 5 or are really a single exponential over 
the full radial range is a question that requires more data 
at R < 5 kpc. 

These present-day patterns in the abundance and spa¬ 
tial structure of the disk must be the consequence of 
the radius-dependent chemical evolution, convolved with 
the subsequent orbit evolution. One possibility is that 
what we are seeing at low [a/Fe] are the (approximate) 
equilibrium points of chemical evolution, where a steady 
state of metal consumption and gas dilution is main¬ 
tained. Most stars in simple chemical-evolution models 
form near the equilibrium [Fe/H], the value of which de¬ 
pends primarily on how much gas is lost to outflows (e.g., 
Binney & Merrifleld 1998; D. Weinberg, et ah, 2016, in 
preparation). Outflows are likely more effective in the 
outer disk than in the inner disk; this scenario would nat¬ 
urally explain the well-deflned radial range spanned by 
each low-[(a/Fe] MAP and the anti-correlation between 
Fe/H] and Rpeak- The radial migration that the low- 
a/Fe] populations have likely experienced (see below) 
will have smoothed the radial profiles such that the ini¬ 
tial radial profile of each low-[(a/Fe] MAP was even more 
sharply peaked around a single radius. 

How do the high-[(a/Fe] MAPs with Rpeak < 5 kpc fit 
into this scenario? The radial profiles of MAPs with 
the same [Fe/H] but different [a/Fe] are strikingly dif¬ 
ferent, especially at low [Fe/H]. This suggests that they 
are not connected by an evolutionary track from high- 
to low-[(a/Fe] at a given [Fe/H]. The similarity in the 
radial profile of high-[(a/Fe] MAPs, combined with the 
narrow range of [a/Fe] spanned at each [Fe/H] along 
the high-[(a/Fe] sequence (Nidever et al. 2014), points 
towards formation and evolution scenarios of all high- 
[a/Fe] MAPs that were similar enough to be structurally 
indistinguishable by now. Because high-[a/Fe] popula¬ 
tions are likely to be the oldest populations in the disk 
(see discussion in Bovy et al. 2012c and Haywood et al. 
2013; Bergemann et al. 2014; Martig et al. 2015), this 
implies that similar physical conditions existed through¬ 
out the disk at early times, a conclusion also reached 
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by Nidever et al. (2014), because of the constant locus 
of the high-[(a/Fe] sequence in the ([Fe/H], [a/Fe]) plane 
throughout the disk. What, if anything, caused this to 
change for the low-[a/Fe] sequence remains to be sorted 
out (e.g., Stinson et al. 2013; Bird et al. 2013). 

We have also discovered that the thickness of low- 
[a/Fe] MAPs is not constant, but instead flares in an 
approximately exponential manner (Figures 12 and 13). 
Such flaring is commonly seen in simulations of the 
outer disk, because of the dynamical heating due to or¬ 
biting satellites and mergers (e.g., Quinn et al. 1993). 
However, we measure the same level of flaring for the 
outer-disk, low-[Fe/H] MAPs as we do for the centrally- 
concentrated, high-[Fe/H] MAPs. It is improbable that 
the latter have been affected much by outer-disk satel¬ 
lite heating. Flaring with an exponential profile and 
at the level that we detect is an important prediction 
of J;^-conserving radial migration—that is, any redistri¬ 
bution of angular momentum that approximately con¬ 
serves the vertical action Jz (e.g., Minchev et al. 2012; 
Solway et al. 2012; Roskar et al. 2013). Therefore, we 
consider the observed flaring of low-[(a/Fe] MAPs as an¬ 
other important indication that radial migration signifi¬ 
cantly affects the distribution of low-[a/Fe] stars. This is 
in addition to the fact that one needs radial migration to 
explain the observed lack of correlation between age and 
metallicity for low-[(a/Fe] stars (e.g., Sellwood & Binney 
2002; Schonrich & Binney 2009) and to explain the re¬ 
versal in the skew of the metallicity distribution function 
when going from the inner Galaxy to the outer Galaxy 
(Hayden et al. 2015). 

We find at > 5cr confidence that the high-[(a/Fe] MAPs 
do not flare in the same manner as the low-[(a/Fe] MAPs. 
Instead, we find that they must have nearly constant 
thicknesses. At face value, this finding implies that little 
radial migration has occurred in the high-[(a/Fe] popula¬ 
tions, because effects such as those from mergers that 
can undo the flaring by mixing in-situ and migrated 
populations (Minchev et al. 2014) do not apply to the 
mono-abundance populations considered here. The fact 
that large-scale migration does happen for the low-[(a/Fe] 
populations entails that whatever causes migration likely 
only strongly affects kinematically-cold populations. Spi¬ 
ral structure whose strength rapidly declines with height 
is an obvious candidate for such a migration mechanism. 
However, determining the exact implications of the con¬ 
stant thickness of high-[(a/Fe] populations for the level 
at which radial migration affects the thicker, high-[a/Fe] 
populations requires chemo-dynamical models that em¬ 
ploy a realistic model for the diffusion of stellar orbits 
due to migration-inducing perturbers. Whatever the case 
may be, the lack of flaring in the high-[(a/Fe] populations 
makes it unlikely that their thickness is due to migration 
or satellite heating. 

While our measurements of the scale heights of MAPs 
are noisier than those of Bovy et al. (2012c), they are 
good enough to confirm the existence and ubiquity of in¬ 
termediate hz components and the smoothness of the 
transition between the traditional “thin” and “thick” 
disks: the vertical structure of the disk cannot be de¬ 
scribed by only two structurally distinct components. 
Overall, the smoothness of the hz transition, the fact 
that the high-[(a/Fe], large-h^ MAPs are centrally con¬ 
centrated, and the level and homogeneity of the flaring 


observed in low-[(a/Fe] MAPs all point toward smooth, 
internal processes dominating the evolution of the MW 
disk. 

Few cosmological simulations exist that can be directly 
compared to our results. A dissection of the radial pro¬ 
file of galactic disks into mono-abundance or mono-age- 
metallicity (a convenient proxy for MAPs in simulations) 
populations has not been attempted in any simulation. 
The analysis of Stinson et al. (2013) comes closest of any 
simulation, because they plot the radial profile of MAPs 
in the MAGIG simulations in their Figure 1 for a few 
MAPs. However, only a narrow radial range around the 
solar circle is displayed, and larger-scale trends were not 
investigated. In agreement with our measurements, they 
find that high-[(a/Fe] MAPs are well described by single 
exponentials in the radial direction. For the low-[Qf/Fe], 
low-metallicity MAPs they find fiat radial profiles, with 
some of them showing a shallow peak. However, their 
low-[(a/Fe], high-metallicity MAPs do not display the 
broken-exponential that we find here, but are instead 
consistent with exponentials. Most other recent simula¬ 
tions, such as those of Minchev et al. (2013), Bird et al. 
(2013), and Martig et al. (2014) only dissect the disk in 
terms of mono-age populations and find centrally-peaked 
radial profiles for all such populations. Using [a/Fe] as 
an age proxy this does not come as a surprise, as discard¬ 
ing the rich [Fe/H] structure in, e.g.. Figure 10, will lead 
to very different radial profiles from those observed for 
MAPs here. The simplified one-dimensional simulations 
of a thin, axisymmetric, gravitationally-unstable disk by 
Forbes et al. (2012) similar to the turbulent disks found 
in cosmological simulations (e.g., Bournaud et al. 2009) 
do lead to donut-like mono-age populations by shutting 
off the gas supply of the inner regions of the disk over 
time. 

The amount of disk flaring for age- or abundance- 
selected populations in these recent simulations varies, 
from little to no flaring for undisturbed disks found by 
Martig et al. (2014) to the strong flaring in the high- 
resolution simulation studied by Bird et al. (2013). That 
the high-[(a/Fe] or old populations have a constant thick¬ 
ness while the low-[a/Fe] or younger populations flare 
significantly has not been seen in any simulation. If 
the flaring is due to migration in these simulations, this 
may indicate that their resolution is not sufficient to dis¬ 
tinguish between the response of kinematically cold and 
warm populations. 

The overall thickness of the disk is likely to be con¬ 
stant due to the mix of centrally-concentrated, thick 
(high-[(a/Fe]) components and more extended, thinner- 
but-fiaring (low-[(a/Fe]) components. This was also re¬ 
cently found in the simulation of Minchev et al. (2015). 

7. CONCLUSIONS 

We summarize our main results as follows: 

• The excellent radial coverage of APOGEE and its 
APOGEE-RG subsample has enabled a detailed investi¬ 
gation of the radial structure of abundance-selected com¬ 
ponents of the disk (MAPs). Any analysis of the spatial 
density distribution of Galactic low-latitude tracers re¬ 
quires explicit accounting for the 3D dependence of in¬ 
terstellar extinction, even if target selections and obser¬ 
vations are in the NIR. We have applied a new likelihood- 
based formalism (Bovy et al. 2016) for determining the 
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radial profiles of stellar tracers in the presence of dust 
extinction. We have been able to obtain good fits to 
the observed star counts, even in regions of very high 
extinction, with simple models for the spatial stellar dis¬ 
tribution. 

• The radial profile of high-[a/Fe] MAPs is consistent 
with a single exponential over the large radial range over 
which they are observed (4 < R/kpc < 14kpc), with no 
sign of a steeper fall-off at large R. Furthermore, all high- 
[a/Fe] MAPs are consistent with having the same expo¬ 
nential S(i?), with a scale length of Hr = 2.2 ± 0.2 kpc. 
This agrees with the results of Bovy et al. (2012c), who 
find a common Hr = 2.01 ± 0.05 kpc for the high-[a/Fe] 
MAPs and those of Nidever et al. (2014), who find that 
the high-[(a/Fe] sequence remains in the same place in 
the ([Fe/H], [a/Fe]) plane throughout the disk. 

• We discovered that the radial surface-density profiles 
S(i?) of low-[(a/Fe] MAPs are complex: they are not a 
single exponential and are not even monotonically de¬ 
creasing outward. Each MAP displays a peak radius 
^peak with an approximately exponential drop-off away 
from i^peak at smaller and larger radii. Thus, the low- 
[a/Fe] stellar disk may be thought of as a sequence of 
narrow, donut-like annuli of increasing i^peak for decreas¬ 
ing [Fe/H]. 

• The peak radius of the low-[a/Fe] MAPs depends 
strongly on metallicity, peaking far inside the solar ra¬ 
dius for the metal-rich low-[(a/Fe] MAPs, and well out¬ 
side the solar circle for the metal-poor low-[a/Fe] MAPs. 
This is consistent with the known radial metallicity gra¬ 
dient. The MAP with solar abundances peaks at the 
solar radius, clearly demonstrating that the Sun is typi¬ 
cal for its Galactic location. 

• The thickness of the high-[Qf/Fe] MAPs is constant with 
R and does not display any flaring. The constraint on 
the inverse flaring scale length of a model with exponen¬ 
tial flaring is strong: R^^^^ = 0.00 ± 0.02kpc~^, when 
combining constraints from multiple MAPs. This argues 
against the local vertical structure of the thick disk com¬ 
ponents being set by outward radial migration. 

• Low-[(a/Fe] MAPs present clear evidence of flaring, 
with an exponential hz{R) profile and a common flaring 
scale length of 8.5 ± 0.7 kpc. This flaring is present both 
for low-[Fe/H], outer-disk MAPs and for high-[Fe/H], 
inner-disk MAPs. 

• We confirm the result of Bovy et al. (2012c), who found 
that all MAPs have a single vertical scale height, with 
a continuous distribution of them from the thinnest to 
the thickest components of the disk. The high preci¬ 
sion abundances from APOGEE (cr[Fe/H] ~ 0.05 dex and 
^[a/Fe] ~ 0.02 dex) renders it unlikely that this smooth 
increase in hz is due to contamination between nearby 
abundance bins. 

The measurements of the vertical profile of MAPs here 
are somewhat noisy, because of a lack of data at inter¬ 
mediate and high latitudes. To make progress on this 
front requires high-latitude data with a known selection 
function. Such data will be provided by the APOGEE-2 
survey (Sobeck et al. 2014). 

We have not attempted to determine updated total 
stellar surface densities E(i?o) or total disk masses asso¬ 
ciated with each MAP here, as was done by Bovy et al. 
(2012a) for the SEGUE G-dwarf sample. Doing so for the 


RG tracers defined using the cuts of Bovy et al. (2014) 
has large uncertainties, because the manner in which gi¬ 
ants trace the underlying population depends on the star 
formation history. Additionally, the severe cuts used to 
define the RG make the sampling of the underlying stel¬ 
lar population highly non-trivial. We plan to determine 
total stellar-population masses for the MAPs in the fu¬ 
ture by using the full APOGEE giant sample and the 
SEGUE G-dwarf sample in combination with the spatial 
densities measured here. 

The current results provide greatly improved con¬ 
straints on the global abundance-spatial distribution 
of stars, which calls for rigorous and global chemical- 
evolution modeling. As age constraints for stars beyond 
the solar neighborhood (^1 kpc) become available for 
extensive samples, generalizing the current analysis to 
include age in addition to abundances in defining MAPs 
will be an obvious step. This opens up the prospect 
of a global map of the Galactic disk in age-abundance- 
position space, even before the release of Gaia data. 
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